Linear density response function within the time-dependent 
exact-exchange approximation 
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We have calculated the frequency-dependent exact exchange (EXX) kernel of time-dependent 
(TD) density functional theory employing our recently proposed computational method based on 
cubic splines. With this kernel we have calculated the linear density response function and obtained 
static polarizabilites, van der Waals coefficients and correlation energies for all spherical spin com- 
pensated atoms up to Argon. Some discrete excitation energies have also been calculated for Be 
and Ne. As might be expected, the results of the TDEXX approximation are close to those of TD 
Hartree-Fock theory. In addition, correlation energies obtained by integrating over the strength of 
the Coulomb interaction turn out to be highly accurate. 
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I. INTRODUCTION 

The present paper is one in a series of papersii^i^ 
reporting on work with the overall aim of finding com- 
putationally efficient but still accurate ways of calculat- 
ing excited-state properties of systems in which excitonic 
effects play an important role. Traditionally, such ef- 
fects have been studied by solving approximate versions 
of the Bethe-Salpeter equation in which the particle- 
hole interaction is taken to be a statically screened 
Coulomb potentiaL* Unfortunately, such ab initio meth- 
ods are computationally very demanding especially in 
low symmetry systems such as nano structures and large 
molecules. During the last decade, time-dependent den- 
sity functional theory^iii^ (TDDFT) has emerged as 
a competing technique due to its computational effi- 
ciency and better scaling with the size of the system. 
In recent years, the use of TDDFT has virtually 'ex- 
ploded' within theoretical chemistry. On the other hand, 
within TDDFT, the limitations are rather caused by 
our rudimentary knowledge of the 'mysterious' exchange- 
correlation (XC) kernel /xc, into which all effects beyond 
the RPA (Random Phase Approximation) are trans- 
ferred. 

In most applications of TDDFT one resorts to an XC 
kernel constructed from some ground-state XC potential 
evaluated at the instantaneous electron density. In this 
way, the non-locality in time (memory effects) leading to 
a frequency-dependent XC kernel, is neglected. These 
are the so-called adiabatic approximations and the sim- 
plest example is the adiabatic local density approxima- 
tion (ALDA) derived from the ground-state local den- 
sity approximation. In fact, any approximate functional 
within ground-state density functional theory (DFT) can 
yield an adiabatic approximation within TDDFT. Al- 
though such approximations have been shown to provide 
good estimates of many physical quantities, there are 
qualitative features which cannot be accounted for. For 
the description of, e.g., excitation energies with multiple- 
particle character the kernel is expected to have a strong 
frequency-dependence^ 



Alternative approaches based on many body perturba- 
tion theory (MBPT) to generate new and more advanced 
kernels have been proposed by some authors j^'^°i^^'^^i^^ 
Up to now, however, the performance of such kernels in 
describing excited- as well as ground-state properties has 
been the subject of little investigation, especially in finite 
systems. 

The simplest approximation derived from MBPT is 
the time-dependent exact- exchange (TDEXX) approxi- 
mation. It can be obtained from a stationary action prin- 
ciple by retaining only the exchange part of the XC part 
of the total action, i.e., terms up to first order in the 
Coulomb interaction. The exact-exchange (EXX) kernel 
/x is frequency-dependent and therefore fundamentally 
differs from the adiabatic approximations. Because the 
EXX kernel can also be derived from our— variational ap- 
proach to the many-body problem, as is done here, it au- 
tomatically posesses conserving properties which means, 
e.g., that the resulting linear density response function 
will obey the /-sum rule. Furthermore, the correspond- 
ing EXX potential of ground-state DFT has already been 
thoroughly studied and shown to share many properties 
of the exact XC potential like, e.g., the correct — l/r 
asymptotic decay and the derivative discontinuity with 
respect to particle number. 

Implementations of the TDEXX approximation has so 
far been limited to the calculation of the total energy 
and the plasmon dispersion relation of the electron gas^ 
the optical absorption spectrum of bulk silicouji^ and, 
recently, van der Waals coefficients and polarizabilities 
for some simple atoms and molecules.— The adiabatic 
TDEXX and the exact TDEXX of a two-electron system 
(which turns out to be frequency-independent) known 
as the PGG (Petersilka, Gossman, and Gross) approxi- 
mation have also been used in the calculation of atomic 
and molecular transition frequencies ;^^i^^ '^^ In the non- 
linear regime the TDEXX approximation has most re- 
cently been applied to the problem of electron dynamics 
in a quantum well.^^ 

In this paper, we calculate the linear density response 
function using the fully frequency-dependent EXX ker- 



nel for all spherical spin compensated atoms up to Argon 
and present results on correlation energies, van der Waals 
coefficients, static polarizabilites and a few discrete ex- 
citation energies for Beryllium and Neon. The correla- 
tion energies, calculated from the HcUman-Feynman the- 
orem applied to the strength of the Coulomb interaction, 
turn out to be very accurate, whereas polarizabilities and 
van der Waals coefficients are similar in quality to the 
rather poor results of time-dependent Hartree-Fock the- 
ory (TDHF). 

The paper is organized as follows. In Sec. [Uwe sketch 
the derivation of the relevant equations and discuss the 
TDEXX approximation in comparison to TDHF. We also 
give a short description of the computational methods. 
In Sec. Illll we present our results and compare them with 
other approximations and exact results. Some attention 
is given to the kernel itself and we provide evidence of 
the /-sum rule being obeyed by studying the large uj- 
behavior of both /^ a-nd the dynamical polarizability. 

As mentioned above, the fact that the TDEXX approx- 
imation obeys the /-sum rule follows from the possibility 
to derive it from the variational and conserving approach 
to MBPT. In the Appendix we show, however, the de- 
tails on how the defining equations lead to the /-sum 
rule. This is rather illuminating and demonstrates the 
necessity for using the correct local exchange potential 
from the LSS equation in the evaluation of the /x kernel 
in order to have the sum-rule fuUfiUed. 

Finally, in Sec. IIV( we draw our conclusions and an- 
nounce a forthcoming publication on spectral properties. 



II. THEORY AND COMPUTATIONAL DETAILS 

A. General formulation 

Within TDDFT the electronic linear density response 
function x is given by 



X = Xs + Xs{v + f,,c)x, 



(1) 



where Xs is the Kohn-Sham (KS) linear density response 
function, v is the Coulomb interaction and /xc is the XC 
kernel defined as the functional derivative of the XC po- 
tential V^r. 



/xc — 



5n 



(2) 



It has been shown in previous publications, ^'A that var- 
ious consistent and, in particular, conserving approxima- 
tions to Wxc and /xc can be obtained from the Klein action 
functional — by choosing physically reasonable approxi- 
mations to the defining $-functionalj^ The stationary 
property of the Klein functional with respect to Green 
functions generated by local potentials leads directly to 
the hnearized Sham Schliiter— (LSS) equation 



j Xs(l, 2)t;xc(2)rf2 = / E,(2, 3)A(3, 2; I)d2d3, 




,/x. 





FIG. 1: Diagrammatic representation of Eq. @ witli the self 
energy in the Hartree-Fock approximation. 



where the self energy Eg is <I>-derivable and expressible 
in only the KS Green function Gs and the Coulomb in- 
teraction. The quantity L is the functional derivative of 
Gs with respect to the KS potential V , 

zA(3,2;l) = ^^^ = G,(3,l)G,(l,2). 

The equation for the corresponding kernel /xc is obtained 
by varying Eq. ^ with respect to V . The result is: 



X.(l,2)/„(2,3),Y.(3.4)(i2<i3 
^A(3.2;l,<i2.3 

A(l,2;4)A(2,3)Gs(3, 1)d2d3 
/"G,(l,2)A(2,3)A(3,l;4)d2d3, (4) 



where A(2,3) = Ss(2,3) - v^^{2)5{2,?,). Due to the 
variational property of the Klein functional and the $- 
derivability of the self energy the kernel /xc obtained from 
Eq. ([4]) will result in a response function (Eq. ([1])) which 
is particle conserving. In the linear regime this means, 
e.g., that it must obey the /-sum rule. 

In this work we are interested in studying the so-called 
TDEXX approximation which is derived at the TDHF 
level of MBPT. The self energy is then given by 



S^(2,3)==z«(2,3)G,(2,3), 
and its variation with respect to V becomes 



'?S;(2,3) 
5V{A) 



-t.(2,3)A(2,3;4). 



(5) 



(6) 



(3) 



The resulting equation for the EXX kernel /x is rep- 
resented diagrammatically in Fig. [T] Notice that the 
obtained /x is often referred to as the kernel of the 
TDOEP (time-dependent optimized effective potential) 
method. It has been derived several times before by 
other people starting with Sharp and Horton^^ in the 
fifties and continuing with Talman and Shadwick^** in 
the sixties. The resulting response function has been 
derived by Rajagopal,^^ Holas, Aravind and Singwi;^ 
Lemmens, Brosens and de Vresej^i^ and more recently 
by Gross^^ to mention a few. Our way of deriving the 



TDHF 



X 



<0*<£>^^^^J^^^^^^^f^^<U>^ 



O 



X.S + Xs/xXs + XsVXs + 




FIG. 2: The linear density response function in TDHF . The first row shows a diagrammatic expansion using the HF Green 
fiinction. The second row shows an expansion in terms of the KS Green function. All diagrams up to first order are drawn and 
seen to be the same as the first order terms of x in TDEXX. 



same approximation has the advantage of demonstrating 
the conserving properties of the approximation. These 
are properties also inherent in the original TDHF ap- 
proximation. Below we will make further comparisons 
between TDHF and TDEXX, illustrating their similari- 
ties and differences. 



B. TDEXX as compared to TDHF 

The quantity to be determined in the TDHF equation 
is the three-point function P = 5G/5Vc-i^ti where T4xt 
is the external potential and G is the HF Green function 
given by the solution of the Dyson equation, 



G = Go + GoS"^ [G]G 



(7) 



where Go is the non-interacting Green function contain- 
ing just the external (nuclear) potential and E^^ — 
Vh + S'^ with Vh being the Hartree potential. By varying 
the Dyson equation with respect to the external poten- 
tial we arrive at the equation for P, usually referred to 
as the TDHF equation 



P(l,3;2) = G(1,2)G(2,3) 



rf4567G(l,4)G(5,3) 



(5SH^(4, 5) 
<5G(6,7) 



P(6,7;2)(8) 



The linear density response function can then be ob- 
tained as x™"^(l,2) = -iP(l,l;2). In Fig. H on the 
first row, x^^^^ is depicted diagrammatically with terms 
up to second order in the explicit dependence on the 
Coulomb interaction. There is, of course, also a depen- 
dence on V through G, which is summed up to infinite 
order in the interaction strength through Eq. ([7]). 

Instead of working with Go as our zeroth order Green 
function we can choose to work with the KS Green func- 
tion Gs- This Green function can also be found through 
a Dyson-like equation: 



G, - Go + Go{Vb_[Gs] + v^[Gs\}Gs 



(9) 



By inverting Eq. ([7]) and Eq. ^ we find again the equa- 
tion for G but with Gg as the zeroth order Green function: 

G^Gs + G,{eHP[G] - V^[Gs] - v^[Gs]}G. (10) 



Iterating to first order gives us 

G(i) = G, + GsimOs] - v^[Gs\}Gs. (11) 

The strictly first order diagrams (in terms of Gs) for 
^TDHF ^^^ j^g.^ i^g identified (see the second row in 
Fig. [2) and we see that they are identical to the first or- 
der terms of x in the TDEXX approximation (see Eq. ^ 
and Fig. [T]). Thus, wc can conclude that, to first order 
in V, TDHF and TDEXX are the same. Notice that, 
in principle, this conclusion is independent of the choice 
of Gs and Wx as long as they are related via Eq. ^. 
But, by using the self-consistent EXX potential the cor- 
responding density is optimized to exactly reproduce the 
HF density up to first order, thus minimizing the con- 
tribution of the approximate higher order terms in x- 
However, using any other reasonable Wx the correspond- 
ing X in TDEXX is still expected to be rather close to the 
response function of TDHF. Indeed, as we shall see later 
(Sec. IIIIFp . using the EXX potential, the LDA potential 
or the exact XC potential leads to very similar response 
functions. The higher order terms of the TDEXX series 
can thus be interpreted as an approximation to the cor- 
responding higher order diagrams of TDHF, where the 
frequency- independent four-point kernel of Eq. ([H]), i.e., 
6Tj^/6G, is simulated by the frequency-dependent two- 
point kernel /x. The trick to approximate the beyond 
first order terms in the series of the TDHF response in 
terms of the zeroth and first order terms such that the 
whole series can be summed as a geometric one has been 
suggested beforci^ The full TDHF series can be written 
order by order as 



^TDHF ^ ^^[1 + ^-1^^+ ^-1^2 + 



= Xo + Xi +XiXo Xi + 



Xo 



1-Xo^Xi 



If Xo = Xs this is just the TDEXX series. 

As shown in the Appendix it is essential to derive the 
potential Wx from the LSS in order to have the /-sum 
rule obeyed. Self-consistency is, however, not necessary. 
Choosing any density there is a local one-body poten- 
tial which, in a non-interacting system, generates that 



density as well as one-electron orbitals with correspond- 
ing eigenvalues. The LSS within EXX (Eq. (O) then 
gives a potential v^ which, together with these orbitals 
and eigenvalues, yields an /x from Eq. Q. In this way, 
both Wx and /x are well defined functionals of the starting 
density. If the density is that given by EXX we have self- 
consistency and the resulting /x is that presented here. If 
we instead start with a better density, much closer to the 
exact one, the corresponding eigenvalue differences will 
be much closer to real particle-hole excitation energies^ 
and we will obtain an /x which will still obey the /-sum 
rule while giving rise to an optical spectrum with a con- 
tinuum starting in almost the correct place - because the 
highest occupied exact DF eigenvalue equals the the neg- 
ative of the ionization potential. "^^ This topic is further 
discussed together with numerical results in Sec. IIIIFI 



C. Numerical method 

The numerical implementation of the TDEXX approx- 
imation starts with a calculation of the self-consistent 
ground-state potential v^ from Eq. Q. This potential 
determines the KS system from which Gg and Xs sue 
calculated and inserted into Eq. (j4]). The kernel is then 
obtained through multiplication by the inverse of Xs both 
from the left and the right. Finally, to obtain the full re- 
sponse function x a R-PA like equation needs to be solved, 
Eq. ([1]). Since /x is frequency dependent all steps after 
the calculation of the potential need to be repeated at 
every frequency. Notice also that, due to the spherical 
symmetry of our studied systems, Eq. ^ separates into 
several decoupled equations, one for each angular mo- 
mentum channel. 

We have chosen a basis set implementation using cu- 
bic splines as radial basis functions. This basis set has 
shown to be ideally suited for solving the LSS equation)^ 
an equation known to be numerically unstable when us- 
ing other methods and basis functions.— i^i^i^i^ We 
have here found that cubic splines also provide an effi- 
cient method for solving the equation for /x. In retro- 
spect this might not be surprising since this equation has 
similarities to the LSS equation. 

A detailed description of the construction of our basis 
set can be found in Ref. Q and a general discussion of B- 
splines in electronic structure calculations can be found 
in Refs. ISalSTI . Here, we will only list the advantages of 
using these basis functions. 1) To start with, they are lo- 
cal functions. This gives us a great amount of flexibility 
in choosing the distribution of splines. Where high accu- 
racy is needed, like in our case close to the nucleus, the 
density of splines can be chosen arbitrarily high without 
loosing accuracy in other regions. 2) Another important 
property of the splines is that there is no risk of instabil- 
ities due to overcompletness because of the strong local- 
ization of the splines. 3) Since every spline only overlaps 
with its three nearest neighbors all matrices will be band 
diagonal, reducing the amount of storage needed and al- 



lowing for the use of efficient diagonalization algorithms. 
4) Once a mesh distribution (e.g. a power law or an expo- 
nential distribution) and a maximum radius is set there is 
only one numerical parameter to vary, i.e., the number N 
of cubic splines. 5) The basis set is complete. The results 
should thus converge to the exact results as A^ — > oo. 6) 
We have shown that the product of two orbitals can be 
re-expanded in the same basis set without increasing the 
number of basis functions. All two-particle functions, like 
for instance response functions, thus become matrices of 
the same order as one particle propagators. 7) A cubic 
spline is composed of cubic polynomials and hence all in- 
tegrals can be solved exactly, either analytically or, and 
actually faster by, using simple Gaussian quadrature. 

The work to calculate the density response function 
within just the RPA (without exchange) is just about as 
extensive as that required to obtain the response function 
from any so called adiabatic - or frequency independent 
- approximation. The calculation of the exchange kernel 
/x of the EXX involves sums over two continua for every 
frequency. An ordinary RPA calculation using N basis 
functions requires for every frequency const * N^ opera- 
tions because of the necessity to invert an N x N matrix. 
Including also an exchange-correlation kernel requires for 
each frequency a double sum over the continuum, i.e., 
const * N'^ operations plus two additional matrix inver- 
sions. This increases the prefactor of the N^ dependence 
on the number N of included splines. The prefector also 
depends heavily on the number of occupied states but, 
for large N, including the exchange kernel does not con- 
stitute a qualitative difference compared to an ordinary 
RPA calculation as far as the calculational effort is con- 
cerned. 

To get an accurate description of both the occupied 
and the first few unoccupied orbitals we used a cubic 
distribution of mesh points in all our calculations. The 
results were converged with '~ 40 splines for He and with 
'~ 60 splines for Ar. 



III. RESULTS 

In this Section we present our results on static po- 
larizabilities, van der Waals coefficients and correlation 
energies for all spin-compensated spherical atoms up to 
Ar and some discrete excitation energies for Be and 
Ne. If not indicated all results are obtained with the 
self-consistent EXX Green function. The convergence 
criterion for the EXX potential was set to \n^ \r) — 
„(fc-i)(r)| < 10-7. 



A. Static polarizabilities 

The static polarizability is defined according to 

a(0) = - / zx{r, r\ uj = 0)z'drdr'. (12) 



TABLE I: Static polarizabilities for some different atoms cal- 
culated in TDEXX, TDHF, RPA and from the KS system, 
(a.u.) 



Atom 


KS 


RPA 


TDEXX 


TDHF 


Lift.'' 


He 


1.487 


1.199 


1.322 


1.322" 


1.38 


Ne 


2.838 


2.234 


2.372 


2.377" 


2.67 


Ar 


16.965 


9.883 


10.737 


10.758" 


11.08 


Be 


81.385 


33.489 


45.648 


45.62* 


37.8 


Mg 


140.26 


60.262 


81.658 


81.60' 


71.53 



"From Ref. M 
''From Ref. |39 
'^From Ref. 



TABLE H; van der Waals coefficients calculated in different 
approximations, (a.u.) 



Atom 


KS 


RPA 


TDEXX 


TDHF 


Lift. 


He 


1.664 


1.171 


1.375 


1.375 


1.458" 


Ne 


7.492 


5.003 


5.506 


5.524" 


6.383" 


Ar 


128 


54.23 


61.88 


61.88" 


64.30" 


Be 


660 


179 


283 


284* 


214^* 


Mg 


1723 


482 


765 


758" 


627'' 



"From Ref. 
''From Ref 
'^From Ref. 
■^From Ref. 
"^From Ref. 



For a system with spherical symmetry only the angular 
momentum channel L = 1 of x contributes. The small 
polarizabilities of the noble gas atoms as compared to 
the alkali earth atoms are due to the large gap in the 
excitation spectrum. On the contrary, Be and Mg have a 
near degeneracy in the HOMO-LUMO gap causing very 
large polarizabilities. 

In Table [J we compare q;(0) calculated in TDEXX with 
a(0) calculated in RPA, TDHF and from the KS response 
function, Xs. Calculating the static polarizability from 
the KS response function provides a reasonable estimate, 
albeit too large compared to the true static polarizabil- 
ity of noble gas atoms. In Be and Mg the error is much 
larger and leads to an overestimation by a factor of two. 
Including interaction effects at the level of RPA reduces 
the KS results for all atoms. However, the RPA polar- 
izabilities are consistently too low as compared to more 
accurate values. Introducing exchange effects at the level 
of TDEXX increases a{0) again leading to an apprecia- 
ble improvement for the noble gas atoms but the error for 
Be and Mg remain roughly the same but with a different 
sign. 

The TDEXX polarizabilities are seen to be very close 
to those of TDHF. This is not surprising as TDDFT 
within the EXX can be considered as a variational so- 
lution to the integral equation of TDHF theory. For a 
two-electron system like He TDEXX is actually identical 
to TDHF. 



B. van der Waals coefficients 

The van der Waals coefficient, or Cg-coefficient, 
tween ion A and B is given by the formula 



be- 



Cr = - 



aA{iuj)aB{ii^)duj, 



(13) 



where aA{iijj) is the dynamic polarizability of ion A cal- 
culated at imaginary frequencies. In Table |TT] the Cg- 
coefficients in the TDEXX approximation are presented 
and compared to the values of the KS system, values cal- 
culated in the RPA and the TDHF approximation, as 
well as accurate values found in the literature. Although 
TDEXX significantly improve over the RPA for He, Ne 
and Ar the Cg coefficients remain too small. The results 
for Mg and Be in TDEXX are too large and give no im- 
provement over RPA in an absolute sense. The TDEXX 
values is again seen to be in good agreement with the full 
TDHF results. This was also noted in Ref. 15 for He, Ne 
and Ar. 



C. Correlation energies 

Using the standard trick based on the Hellman- 
Feynman theorem to integrate the interaction energy 
with respect to the strength of the Coulomb interaction 
the correlation energy becomes 



Er_ 



dX 

Jo 



dw 

2^ 



Tv{v[x\iu)-Xs{iu;)]}. (14) 



In Eq. (IMI) we have use the short hand notation Tr fg = 
J drdr' f{r,r')g{r' ,r) for any two-point functions / and 
g, and defined the response function 



X 



Xs 



1 - (A« + /4)x. 



where f^^ is the XC kernel of a system of electrons inter- 
acting through the rescaled Coulomb potential Iv. From 
this expression the simplest approximation is obtained 



TABLE HI: Correlation energies calculated within different 
approximations and compared to accurate CI calculations. 
The correlation energy is here defined as the total energy mi- 
nus the HF energy. In the last column also the Hartree-Fock 
total energies are tabulated, (a.u.) 



Atom TDEXX 


RPA 


MP2" 


CI' 


HF'' 


He 


0.044 


0.083 


0.047 


0.0420 


2.8617 


Ne 


0.389 


0.596 


0.480 


0.3905 


128.5471 


Ar 


0.721 


1.091 


0.844 


0.7225 


526.8175 


Be 


0.102 


0.181 


0.124 


0.0943 


14.5730 


Mg 


0.445 


0.681 


0.514 


0.4383 


199.6146 


"From Ref 


131|- 










''From Ref 


n. 











RPA 




Er = E. 



FIG. 3: Diagrams contained in the correlation energy func- 
tional, Eq. (fT4)l . with /xc — /x- 



by setting /xc = 0, leading to the formula for the RPA 
correlation energy. A fully self-consistent calculation of 
this approximation was performed in Ref. 4 for all spin- 
compensated spherical atoms up to Ar. With /xc = /x 
a cancellation between Hartree and exchange terms is 
expected to occur, improving the largely overestimated 
RPA values. We have here, for the first time, performed 
such a calculation. The results are presented in Ta- 
ble mil and compared with those in the RPA and the 
MP2 approximation'^^ as well as results from accurate CI 
calculations.''^ As expected, the TDEXX results are very 
accurate. A diagrammatic analysis (see Fig. ([3])) shows 
that with this kernel the correlation energy will, apart 
from the RPA or bubble series of diagrams, also contain 
the important second order exchange diagram (included 
in the MP2 approximation) as well as an infinite series 
of terms simulating the higher order exchange diagrams. 



D. /-sum rule and EXX kernel 

In the Appendix we prove that the /-sum rule, a con- 
sequence of particle conservation, is valid in the TDEXX 
approximation by showing that the coefficient of the l/uj"^ 
term in the large oj expansion of x is the same as in the 
expansion of Xs- Thus, by expressing the dynamical po- 
larizability, Eq. (fT^ . as 



a(*w) = 51 ■ 



/, 



(15) 



where q = {k,fj,) is a particle-hole index, 
citation energy and fq is the corresponding oscillator 
strength, the sum over all oscillator strengths must equal 
the number of particles 



E/. = ^- 



As a test of our numerical accuracy this result was 
checked. We multiplied a{iuj) by w^ and studied the 
large w- values of this function. Indeed, within high accu- 
racy (^ lO"**), the results converged to N for all atoms. 

In Fig.lHwe have plotted Lo'^a{iui) for He. The kernel /x 
was calculated using either the Green function Gs and the 
exchange potential Vx of the EXX or the same quantities 
within LDA. The /-sum rule is not expected to be obeyed 
in the latter case. And this is, indeed, what is observed. 



1.96 



1.92 




a; [a.u.] 



FIG. 4: The main figure shows LO'^a{iLo) for He calculated at 
different Green functions. The inset shows /x' with L = 1 and 
q = 2s ~* 2p for Be. The large w-behavior in both plots clearly 
indicates that the /-sum rule is obeyed in our calculations. 



It should be noticed, however, that the violation is minor 
(- 0.8%). 

The quantity XsfxXs must decay as 1/w'* (see the Ap- 
pendix). As a consequence, /x cannot diverge (as it does, 
e. g., in a model system, see Ref. |40) and must approach 
a constant as w — > oo. Another way to check the /-sum 
rule and to test our calculations is thus to study the large 
w-behavior of /x. To do this we first notice that the ker- 
nel only appears in the form of matrix elements of the 
Fg-functions: 



fr'i^) =jFqir)Ucir,r',^)Fq,{r')d 



rdr' 



where Fq is a KS excitation function, i. e., a product 
of the occupied KS orbital ipk and the unoccupied KS 
orbital (/j^. Now, since the excitation functions alone in- 
tegrates to zero the kernel is unique only up to the addi- 
tion of two arbitrary functions, 5i(w,r) and g2{u},r'). 
The quantity f^'' is unique though and in Fig. [5] we 
have plotted this quantity for Be at imaginary frequencies 
and different L. The frequency dependence is approxi- 
mately constant for low frequencies but becomes more 
pronounced for higher frequencies. This justifies the use 
of the adiabatic approximation for low energies. At large 
u; every matrix elements /^"^ (iuj) is seen to approach a 
constant. This is again demonstrated for one matrix el- 
ement in the inset of Fig. [H 



E. Excitation energies 

The poles of the exact linear density response func- 
tion correspond to the particle conserving excitation en- 
ergies of the system. Thus, any approximate x yields an 
approximate set of excitation energies. Xs contains the 
excitation energies of the KS system and in many cases 



they give a good approximation to the discrete part of 
the spectrum. Nevertheless, Xs is a non-interacting re- 
sponse function and as such lacks many features of the 
full many body x- 

With /xc = 0, Eq. ^ can be rewritten as an eigenvalue 
problem where the eigenvalues correspond to the poles of 
X- Due to the frequency dependence of the XC kernel of 
the EXX this technique to obtain the excitation ener- 
gies cannot be applied in a straight-forward way. In the 
present work we have chosen to work with localized basis 
functions which causes Imx as a function of frequency 
to consist of a series of sharp delta functions. And this 
behavior is independent of whether the frequency lies in 
the continuum or in the discrete part of the spectrum. In 
order to obtain something which can be plotted we arbi- 
trarily add a small positve imaginay part to the real fre- 
quency and obtain the discrete excitation energies from 
the positions of the resulting huge peaks in Imx- 

In Table HV] we present the first few discrete excitation 
energies of Be and Ne and compare them to experimental 
values and to the ones obtained in RPA, TDHF and the 
KS system (here, meaning the differences between the KS 
one-electron eigenvalues). In the case of Be the KS val- 
ues are too low. The RPA improves over the KS results 
but proceeding to the TDEXX makes them, somewhat 
worse. It is interesting to observe how close the results 
of the TDEXX are to those of the TDHF. In the case of 
Ne we have observed a qualitatively different behavior. 
The KS eigenvalue differences are larger than the exper- 
imental values. As in the case of Be, the RPA tends to 
increase the excitation energies yielding, for Ne, an even 
larger discrepancy. The inclusion of exchange effects does 
not lower the RPA results, as in the Be case, but rather 
increase them. Although TDEXX and TDHF push the 
RPA excitation energies in the same direction the actual 
values are not as close as for Be. 



0.04 




FIG. 5: The quantity /^'' for Be is plotted for different L at 
imaginary frequencies. There are two curves with L = 1. In 
the lower, q corresponds to the 2s — > 2p transition and in the 
upper to the Is — > 2p transition. 



The oscillator strengths for the discrete excitation en- 
ergies can be extracted from the height of the peaks in 
the optical spectrum. As an example, we obtain 1.379 
for the lowest transition of Be. This can be compared to 
the TDHF value of 1.37841 



F. Sensitivity to the ground state KS 
approximation 

As noticed by otherSfl^, it is crucial to have good KS 
transitions as a starting point when calculating the dis- 
crete excitation energies. With the exact KS potential^ 
for Ne the 2p -^ 3s transition in the KS system be- 
comes 0.612, an already very good approximation to the 
true value (0.619). The transition frequencies we have 
obtained are thus expected to be improved by using a 
ground-state potential better than that of EXX. When 
we improve the ground-state density and the correspond- 
ing orbitals and eigenvalues we still calulate v^ from the 
LSS in order not to violate the sum-rule as discussed in 
Sec. imi 

We have used the exact densities of Umrigar et. al4^ 
and found the static polarizabilities 1.35 and 40.5 for He 
and Be, respectively. There is, thus, a significant im- 
provement, particularly in the case of Be. This shows 
that the largest error is actually caused by a poor de- 
scription of the ground state within the EXX. With the 
very accurate densities by Umrigar et. al. the transition 
2s -^ 2p in Be becomes 0.182 and the 2p -^ 3s transition 
in Ne becomes 0.631. The transition frequencies are thus 
also largely improved by using an accurate KS ground 
state. 

By also using the accurate XC potentials by Umrigar 
et. al. in Eq. (|4]) for generating /x we violate the /- 
sum rule but, in this way, as discussed in Sec. IIIBl we 
obtain a response function which to first order in the 
Coulomb interaction is identical to that of TDHF. We, 
therefore, expect the results to revert to our original re- 
sults from TDEXX which are close to those of TDHF. 
Indeed, we now obtain the static polarizabilities 1.323 
and 45.76 for He and Be, compared to our previous re- 
sults 1.322 and 45.65 respectively. Using instead the LDA 
potential for generating /^ we obtain 1.343 and 45.23 for 
the same quantities, again demonstrating the closeness 
to the TDHF. 



IV. CONCLUSIONS AND OUTLOOK 

In this paper we have calculated the linear density re- 
sponse function of the TDEXX approximation for all 
spin-compensated spherical atoms up to Ar. For the 
properties studied in this work, i. e., static polarizabil- 
ities, van der Waals coefficients and the low-lying exci- 
tation energies, the results show that TDEXX is a good 
approximation to TDHF. 
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TABLE IV; The first few discrete excitation energies for Be 
and Ne in TDEXX compared to experimental, TDHF, RPA 
and KS transitions. 



Transition 


KS 


RPA 


TDEXX 


TDHF" 


Exp.* 


Be 

2s--*2p 


0.1312 


0.2032 


0.1764 


0.1764 


0.1940 


2s^3p 


0.2412 


0.2547 


0.2470 


0.2471 


0.2742 


2s^4p 


0.2731 


0.2777 


0.2749 


0.2750 


0.3063 


2s-»5p 


0.2868 


0.2889 


0.2877 


0.2878 


0.3195 


Ne 












2p^3s 


0.6585 


0.6675 


0.6803 


0.6739 


0.6190 


2p-^4s 


0.7793 


0.7812 


0.7827 


0.7818 


0.7268 


2p^5s 


0.8134 


0.8141 


0.8147 


0.8139 


0.7593 














"From Refs. \MM- 




''Adopted from Refs. H^EI. 





TDEXX only takes into account exchange effects in the 
response of the system to external perturbations. For no- 
ble gas atoms the static polarizabilities and van der Waals 
coefficients still turn out to be in rather good agreement 
with the experimental values. The relative error is about 
5% for He, 13% for Ne and 3% for Ar. On the other hand 
for alkali earth atoms the results are not as satisfactory, 
the relative error being around 26% for Be and 18% for 
Mg. For these systems it is thus necessary to include 
correlation effects in order to get an accurate description 
of the above properties. The low-lying excitation ener- 
gies are accurate to within 10% for both Be and Ne. We 
show here, however, that significantly better excitation 
energies as well as static polarizabilities can be obtained 
by using a more accurate exchange correlation potential 
for the ground state. 

We have also calculated the total energies from the 
Hellman Feynman theorem applied to the strength of the 
Coulomb interaction. The results are in excellent agree- 
ment with accurate CI calculations. This is probably due 
to the fact that the fluctuation-dissipation formula at the 
level of TDEXX accounts for several correlation diagrams 
as, for instance, the important second order exchange di- 
agram. 

Finally, we have examined the behavior of kernel /x 
along the imaginary frequency axis and found it to be 
both weakly and slowly dependent on u. It is, how- 
ever, not without structure indicating a more complex 
behavior on the real axis as compared to that resulting 
from simple poles. Indeed, we have seen that the ker- 
nel has both single and double poles on the real axis. 
The weak frequency dependence at the imaginary axis 
suggests that the adiabatic, i.e., frequency independent 
approximations to the EXX approximation, might not be 
such a bad idea at least not as far as total energies are 
concerned. The experience from the electron gas, how- 
ever, strongly contradicts this conjecturei^ 

In most cases of approximations within DFT and 
TDDFT it is very difficult to see through what kind 
of physical processes are actually incorporated into that 
particular approximation. In our case, basing our ap- 



proximations on <I>-derivable theories within MBPT, we 
can say with confidence that a description of double exci- 
tations is way above the EXX. Such a description would 
require an /xc based at least on the time-dependent GW 
approximation which, by the way, probably would yield 
much better van der Waals coefficients. But again we 
refer this discussion to a future publication. 
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APPENDIX A: /-SUM RULE 

In Ref. 13 we demonstrated that every /xc of TDDFT 
obtained from the variational formulation of that work 
obeys particle conservation, which amounts to the /-sum 
rule in the linear limit. 

In this Appendix we will demonstrate how the /-sum 
rule explicitly comes out of the construction of the ker- 
nel /x of the EXX approximation within TDDFT. Among 
other virtues, this detailed derivation demonstrates the 
crucial importance of using the LSS equation when con- 
structing the kernel /x in order to have the sum rule ful- 
filled. From Eq. ([T]) we see that the total density response 
function x can be written as 

X = [1 - {v + u)xsr'^ Xs 

Since the ground-state KS theory gives the correct den- 
sity it follows that Xs obeys the /-sum rule. Given the 
Lehmann representations for the response functions x 
and Xs this ensures that the coefficient of the l/w^ term 
of the large uj expansion of Xs has the correct value. If 
the full response function x is also to obey the /-sum rule 
there must, obviously, be no contribution to the l/w^ co- 
efficient from the denominator which thus must tend to 
unity at large co. A sufficient condition for this is that 
fxXs vanishes as 1/uj'^ in this limit. From Fig. [T] we see 
that the quantity which is actually calculated from the di- 
agrams is XsfxXs which thus should decay as l/w* given 
the l/w^ of dependence of Xs in this limit. In order to 
see that this is indeed the case we explicitly examine the 
large frequency behavior of the different contributions 
to XsfxXs ■ As seen from Fig. [T] this quantity has con- 
tributions from five Feynman diagrams of which those 
proportional to the exchange potential v^ are most eas- 
ily combined with one each of the two self-energy dia- 
grams. The contribution from one of the diagrams with 



self-energy insertions is 

Si{r, r';u;) ^2i J ^Gs{ri,r;cj') Gs{r',r2;cj') 
xGs{r,r';uj + uj')A{ri,r2) 

where A = v{r,r')J2k'"'k'Pk{r)(pl{r') + v^{r)S{r - r') 
and Hk is 1 for occupied states and zero otherwise. Due to 
time-reversal symmetry the non-interacting Green func- 
tions are symmetric in there spatial arguments. It is then 
easily seen that the contributions from the remaining two 
diagrams with self energy insertions is obtained by adding 
to 5i(u;) above the result 5i(— w). Consequently, the 
sum of all diagrams with self energy insertions is an even 
function of uj. Carrying out the frequency integrations 
we obtain 



5i(r,r';z) 



J2 {k2\A\h)^k, {rM, (r')^fe. {rM, (r') 

kik2k3 



nki - nk3 



nkt - "-fe2 



^k^ ^k2 K ^ ' ^k\ ^k^ ^ I £ki ^k2 

Here, the functions </?fe(r) are the KS orbitals with eigen- 
value Ek and wc have switched from time-ordered quan- 
tities to retarded ones by everywhere replacing lu — i6 hy 
Lo -\- i6. The quantity z is a complex frequency in the 
upper half plane. 

When z ^ oo the contribution appears to behave as 
1/z but wc remind the reader that we should add a sim- 
ilar expression with z replaced by ~z. Then, the leading 
order becomes A/z'^ with the coefficient 



A^A }^ {h\A\h)^k,{rM,{r')ipk2{rM,{r') 

kik2k3 



£k3 - Efca 



Sks ^ £k2 



where we have multiplied by two to account for the re- 
maining two self-energy like diagrams {Si{—io)). We can 
regroup the terms and write A — Ai + A2 where 



Ai =4 



kik2k3 



{k2\A\h)ipk, {r)ipl {r')ipk2 {rM, {r') 



A2-4 ^ {k2\A\hWkArMA^>k2{rM3{r') 



kik2k3 



x(nfc3 -rifcj 



^k\ ^k2 



^k3 ^k2 

Let us now define the one-particle density matrix by 
n{r,r') = 2^nk(pk{r)(p*k{r') 

k 

and use the completeness of the KS orbitals. We obtain 

Ai — 6{r — r') / (i'^r3u(r — r3)|n(r3,r)|^ 

-w(r-r')|n(r,r')P (AI) 



and see that the diagrams containing v^ do not contribute 
to leading order in the large frequency limit. 

In order to manipulate the A2 coefficient we use the 
fact that the KS orbitals obey the KS equation 



__V2 + V{r) } ipk{r) = £fe(^fc(r-) 



where V{r) is the full KS potential. Because of the dif- 
ference Eki — Sk2 , the terms involving the potential will 
vanish and we obtain 



A2^2 J2 {k2\A\k,)^l{r')^l{r 



kik2k3 



,^ nfc3 - nk2 

^k3 ^k2 



X {fk^ (r)V2(^fe, (r) - ipk2 (r)VVfei (r)} 
Now, using 

(flV^(f2 - (^2V^(^i ^ 2V(</?iV(^2) - V^((/?i</?2) 

and the completeness we can write A2 as 

A2 ^ 2V [d{r - r')V/(r)] - 2V^ [S{r - r')f{r)] 
where 



fir) = J2{k2\A\k,)^k2irW*k3ir 



k2k3 



nks ~ nk2 

^ks ^k2 



This is because, by symmetry. 



V.f{r) = 2j2{k2\A\k3Wlir)V^k2ir) 



k2k3 



nk3 - nk2 

£k3 - £k2 



But, /(r) = is the LSS equation defining Ux and we 
have shown that A2 — 0. 

Let us now study the remaining fifth diagram - the 
vertex diagram - in the high frequency limit. The contri- 
bution Rv from this diagram is 

xGs{r,ri;uJi + w)Gs(ri, r'; ^2 + uj) 
xv{ri,r2)Gs{r', r2;uj2)Gs{r2, r; wi) 

Carrying out the frequency integrals and converting to 
the retarded propagator in the upper half plane (uj —* z 
with Imz > 0) gives 

Rv{r,r';u) = -2 ^ "^ Vk^ir)ip*k^{r)lpl,^{r')lf>ki,{r') 

kik2 k'-^k'^ 

n ,M 1;/; \ ('^'=1 -nk2){nk[ - ny^ 

x{kikU\v\k\k2)-, T-, — r 

^ " ' ' '{z + ek2-£k,){z + ek'^-ekO 

Here, the standard Coulomb integral is given by 

{k,k'2\v\k[k2) = j d'rd'r'^l (r)^*, (r')«(r - r') 

xipk'^{r)ipk2{r') 
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In the high frequency hmit, to leading order, this becomes 
Bjuj^ where the coefficient B is given by 

i? = -2 E E ^^1 (^)^fe. (^)^fci ^"-^K (^') 

Using again the completeness and the definition of the 
density matrix n(r^r'^ we obtain 

B = -5{r - r') I Sr^v{r - r3)|n(r3, r)p 
+v{r — r')\n{r, r')\^ 



This is just the coefficient A, Eq. (jAip . with the oppo- 
site sign. Consequently, the 1/z^ contribution from the 
vertex diagram exactly cancels the same term from the 
self-energy contributions meaning that XsfyiXs decays as 
1/z^ at large frequencies. And this provides an explicit 
proof of the /-sum rule in the EXX approximation. 

Finally we note that this result means that the ex- 
change kernel should have a very weak dependence on 
frequency at large frequencies. 
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